% Brief Introduction:
% Spline1DBasis Computes 1-D polynomial spline basis which is specified by
% the structure SplineSpace. This function is simplified from the function
% written by Paul L. Fackler & Mario J. Miranda.
% Usage:
%   [B]=Spline1dBasis(SplineSpace,x,order);
% Inputs:
%   breaks      :    Vector of break points, accendingly sorted
%   k           :    Polynomial order of the spline
%   x           :    Vector of the evaluation points 
%   order       :    The order of differentiation (default: 0)
%                    if a vector, SPLIBAS returns a cell array 
%                    otherwise it returns a matrix
% OUTPUTS
%   B           :    a kxn basis matrix
%
% Note: The number of basis functions is n=size(breaks,1)+k-1
% Author:       Xing Guo, University of Michigan, xingguo@umich.edu
% Version:      Nov 28, 2017               
% USES: Spline1dBasisDoP



function [B]=Spline1dBasis(breaks,k,x,order,SelfLookup)

% For Test Purpose
  
% breaks   =   (0:1:10)';
% k        =   3;
% x        =   (0:0.5:10)';
% order    =   1;

% Parse the Spline Space Structure
% breaks      =   SplineSpace.breaks;
% k           =   SplineSpace.k;

p           =   size(breaks,1);
m           =   size(x,1); 
if nargin<3
    error('At least 3 input should be provided');
elseif nargin==3
    order   =   0;
else
    if isempty(order)
        order   =   0;
    end
end
MinOrder    =   min(order);
MaxOrder    =   max(order);

% A few Checks
if k<0
    error(['Incorrect value for spline order (k): ' num2str(k)]);
end
if min(size(breaks))>1
    error('''breaks'' must be a vector');
end
if MaxOrder>k
    error('Order of differentiation must not be larger than k');
end
if size(x,2)>1
    error('x must be a column vector')
end
% Augment the breakpoint sequence 
n           =   p+k-1;  
a           =   breaks(1); 
b           =   breaks(end);
augbreaks   =   [a(ones(k-MinOrder,1));...
                 breaks(:);...
                 b(ones(k-MinOrder,1))];

%   The following lines determine the maximum index of 
%   the breakpoints that are less than or equal to x,
%   (if x=b use the index of the next to last breakpoint).
if nargin<5
    ind     =   lookup(augbreaks,x,3);
else
    ind     =   SelfLookup(x);
end
% Recursively determine the values of a k-order basis matrix.
% This is placed in an (m x k+1-order) matrix
bas     =   zeros(m,k-MinOrder+1);
bas(:,1)=   ones(m,1);
B       =   cell(length(order),1);

if MaxOrder>0
    D=Spline1dBasisDoP(breaks,k,max(order)); % Derivative op
end 
if MinOrder<0
    I=Spline1dBasisDoP(breaks,k,MinOrder);   % Integral op
end     
for j=1:k-MinOrder
    for jj=j:-1:1
        b0  =   augbreaks(ind+jj-j);          
        b1  =   augbreaks(ind+jj);
        temp=   bas(:,jj)./(b1-b0);
        bas(:,jj+1) =   (x-b0).*temp+bas(:,jj+1);
        bas(:,jj)   =   (b1-x).*temp;
    end
% as now contains the order j spline basis
    ii  =   find((k-j)==order);
    if ~isempty(ii)
        ii  =	ii(1);
      % Put values in appropriate columns of a sparse matrix
        r   =   (1:m)'; r   =   r(:,ones(k-order(ii)+1,1));
        c   =   (order(ii)-k:0)-(order(ii)-MinOrder); 
        c   =   c(ones(m,1),:)+ind(:,ones(k-order(ii)+1,1));
        B{ii}=  sparse(r,c,bas(:,1:k-order(ii)+1),m,n-order(ii));
      % If needed compute derivative or anti-derivative operator
        if order(ii)>0
            B{ii}   =   B{ii}*D{order(ii)};
        elseif order(ii)<0
            B{ii}   =   B{ii}*I{-order(ii)};
        end
    end
end

if MaxOrder==k
    ii  =   find(k==order);
    if ~isempty(ii)
        ii  =   ii(1);
        B{ii}   =   sparse((1:m)',ind-(k-MinOrder),ones(m,1),m,n-MaxOrder);
        B{ii}   =   B{ii}*D{MaxOrder};
    end
end
if length(order)==1, B=B{1}; end